function [output cache] = functions_nsym_polar_geometry(fname, fobj, model, eta, zeta, data, params, cache)

% Set of functions that defines geometrical objects on the quotient manifold 
% St(n,p)\times Cone(p)\times St(m,p)/O(p) where O_p is the orthogonal group
% St(n,p): Stiefel manifold
% Cone(p): Positive definite matrices of size p\times p
% Usages:
%
% Metric:                metric
% Projection on H:       project
% Retraction:            retraction
% Cost function:         f
% Riemannian Gradient:   grad_f
% Riemannian Hessian:    hessian
%
% model is the current iterate
% eta is a search direction
% params is a cell array that contains {I,J,trueDists,EIJ}
% cache is a cell array to reduce redundant computations
%
% Authors:
% Bamdev Mishra and Gilles Meyer
% {b.mishra,g.meyer}@ulg.ac.be



switch lower(fname),
  
    case 'metric', % Riemannian metric
          
        output = trace((eta.U'*zeta.U)) + trace((eta.V'*zeta.V)) + trace((model.B\eta.B)*(model.B\zeta.B));

    case 'project', % Projection onto horizontal space
      
        SS = model.B*model.B;
        AS = model.B*(skew(model.U'*eta.U) + skew(model.V'*eta.V) - 2*skew(model.B\eta.B))*model.B;
        omega = sylvester1(SS,AS);

        projDir.U = eta.U - model.U*omega - model.U*symm(model.U'*eta.U);
        projDir.B = symm(eta.B) - (model.B*omega - omega*model.B);
        projDir.V = eta.V - model.V*omega - model.V*symm(model.V'*eta.V);
            
        output = projDir;
    
    case 'retraction', % Chosen retraction

        if model.p == 1
            model.U = model.U + zeta * eta.U;
            model.V = model.V + zeta * eta.V;
            model.U = model.U/sqrt(model.U'*model.U);
            model.V = model.V/sqrt(model.V'*model.V);            
            model.B = model.B*exp(zeta*eta.B/model.B);
            
        else
            L = chol(model.B);
            model.B = L'*expm(L'\(zeta*eta.B)/L)*L;
            model.U = uf(model.U + zeta * eta.U);
            model.V = uf(model.V + zeta * eta.V);
        
        end
            
        output = model;
        
    case 'f', % Cost function
       
        [output,cache] = feval(fobj, 'f', model, [], data, params, cache);              
        
    case 'grad_f', % Gradient
      
        % No further projection step is required for the gradient
        [output,cache] = feval(fobj, 'grad_f', model, [], data, params, cache);
        
    case 'hessian', % Hessian
      
        % Hessian in the total space
        [hessian,cache] = feval(fobj, 'hessian', model, eta, data, params, cache);
        
        % Hessian on the quotient manifold
        [output,cache] = functions_nsym_polar_geometry('project', fobj, model, hessian, [], data, params, cache);

    case 'zero_dir', % Compute a zero search direction

      sdir.U = zeros(size(model.U));
      sdir.B = zeros(size(model.B));
      sdir.V = zeros(size(model.V));
              
      output = sdir;    

    case 'scale_dir', % Scalar scaling of a search direction

      eta.U = zeta * eta.U;
      eta.B = zeta * eta.B;
      eta.V = zeta * eta.V;
      
      output = eta;
      
    case 'add_dir', % Add two search directions

      model.U = model.U + zeta * eta.U;
      model.B = model.B + zeta * eta.B;
      model.V = model.V + zeta * eta.V;
      
      output = model;
 
end

